options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

hierarchical model

class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)
n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
  scale_shape_manual(values=1:k)

qplot(x,y,col=c)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

estimate as no class

ex8-0.stan

y~N(b00+b10*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-0.stan') 
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)
## Initial log joint probability = -456.673 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23       -165.36    0.00977385    0.00137419           1           1       24    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__    -165.36
##    b0        15.04
##    b1         1.57
##    s         16.56
##    m0[1]     19.24
##    m0[2]     30.59
##    m0[3]     38.28
##    m0[4]     25.01
##    m0[5]     15.61
##    m0[6]     18.88
##    m0[7]     36.08
##    m0[8]     38.06
##    m0[9]     31.31
##    m0[10]    39.85
##    m0[11]    46.01
##    m0[12]    33.20
##    m0[13]    42.04
##    m0[14]    31.87
##    m0[15]    25.03
##    m0[16]    25.33
##    m0[17]    43.10
##    m0[18]    32.24
##    m0[19]    41.91
##    m0[20]    35.20
##    m0[21]    30.50
##    m0[22]    21.31
##    m0[23]    41.22
##    m0[24]    45.78
##    m0[25]    45.40
##    m0[26]    43.43
##    m0[27]    21.12
##    m0[28]    18.46
##    m0[29]    23.21
##    m0[30]    25.39
##    m0[31]    39.09
##    m0[32]    40.21
##    m0[33]    33.91
##    m0[34]    23.15
##    m0[35]    41.73
##    m0[36]    24.18
##    m0[37]    35.80
##    m0[38]    43.58
##    m0[39]    28.55
##    m0[40]    31.78
##    m0[41]    15.39
##    m0[42]    42.26
##    m0[43]    22.36
##    m0[44]    25.72
##    m0[45]    37.82
##    m0[46]    40.62
##    m0[47]    40.85
##    m0[48]    26.78
##    m0[49]    34.30
##    m0[50]    22.82
##    y0[1]     21.34
##    y0[2]     56.24
##    y0[3]     16.27
##    y0[4]     -9.63
##    y0[5]     -6.93
##    y0[6]     38.76
##    y0[7]     37.65
##    y0[8]     43.24
##    y0[9]      3.28
##    y0[10]    41.76
##    y0[11]    48.98
##    y0[12]     1.68
##    y0[13]    30.40
##    y0[14]    42.67
##    y0[15]    -5.45
##    y0[16]    53.63
##    y0[17]    42.79
##    y0[18]    49.27
##    y0[19]    52.37
##    y0[20]    44.08
##    y0[21]    60.54
##    y0[22]    30.18
##    y0[23]    48.02
##    y0[24]    45.54
##    y0[25]    28.12
##    y0[26]    63.27
##    y0[27]    23.60
##    y0[28]    45.35
##    y0[29]    34.86
##    y0[30]    42.53
##    y0[31]    32.93
##    y0[32]    30.24
##    y0[33]    44.82
##    y0[34]    16.25
##    y0[35]    53.02
##    y0[36]    19.29
##    y0[37]    25.12
##    y0[38]    31.93
##    y0[39]    43.55
##    y0[40]    31.07
##    y0[41]    16.88
##    y0[42]    44.27
##    y0[43]     8.55
##    y0[44]    25.93
##    y0[45]    38.93
##    y0[46]    44.69
##    y0[47]    62.46
##    y0[48]    59.17
##    y0[49]    15.37
##    y0[50]    49.14
##    m1[1]     15.39
##    m1[2]     18.79
##    m1[3]     22.20
##    m1[4]     25.60
##    m1[5]     29.00
##    m1[6]     32.40
##    m1[7]     35.80
##    m1[8]     39.20
##    m1[9]     42.61
##    m1[10]    46.01
##    y1[1]     -0.23
##    y1[2]     22.26
##    y1[3]     22.35
##    y1[4]     32.95
##    y1[5]     12.85
##    y1[6]     69.10
##    y1[7]     54.09
##    y1[8]     31.67
##    y1[9]     44.55
##    y1[10]    49.48
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -164.17 -163.85  1.35  1.08 -166.81 -162.74 1.01      439      619
##    b0       14.52   14.54  5.84  5.78    4.43   23.33 1.02      216      256
##    b1        1.60    1.61  0.47  0.47    0.86    2.40 1.02      260      396
##    s        17.38   17.24  1.87  1.82   14.57   20.61 1.00     1984     1300
##    m0[1]    18.82   18.87  4.74  4.67   10.62   25.95 1.02      218      247
##    m0[2]    30.42   30.38  2.54  2.46   26.08   34.44 1.01      403      656
##    m0[3]    38.27   38.31  2.91  2.73   33.36   42.97 1.00     1402     1357
##    m0[4]    24.71   24.80  3.39  3.25   18.73   29.89 1.01      239      279
##    m0[5]    15.11   15.14  5.69  5.65    5.27   23.71 1.02      216      248
##    m0[6]    18.45   18.53  4.83  4.77   10.11   25.74 1.02      217      248
##    m0[7]    36.03   36.02  2.61  2.50   31.67   40.34 1.00     1459     1316
##    m0[8]    38.05   38.07  2.87  2.68   33.17   42.76 1.00     1423     1325
##    m0[9]    31.15   31.14  2.49  2.44   26.95   35.13 1.00      467      785
##    m0[10]   39.88   40.01  3.19  3.00   34.52   45.13 1.00     1146     1403
##    m0[11]   46.17   46.30  4.58  4.40   38.39   53.70 1.01      578     1338
##    m0[12]   33.09   33.10  2.44  2.38   29.01   37.05 1.00      749     1049
##    m0[13]   42.12   42.26  3.64  3.45   35.96   47.97 1.01      857     1363
##    m0[14]   31.73   31.75  2.46  2.41   27.59   35.72 1.00      534      832
##    m0[15]   24.74   24.82  3.38  3.24   18.76   29.92 1.01      239      288
##    m0[16]   25.04   25.11  3.32  3.18   19.19   30.18 1.01      242      289
##    m0[17]   43.21   43.34  3.88  3.65   36.67   49.47 1.01      752     1303
##    m0[18]   32.11   32.12  2.45  2.38   27.97   36.10 1.00      580      865
##    m0[19]   41.99   42.13  3.61  3.41   35.86   47.80 1.00      874     1363
##    m0[20]   35.13   35.09  2.53  2.41   30.87   39.26 1.00     1312     1306
##    m0[21]   30.33   30.29  2.54  2.47   25.97   34.39 1.01      396      656
##    m0[22]   20.93   20.98  4.22  4.16   13.63   27.29 1.02      220      243
##    m0[23]   41.28   41.39  3.46  3.29   35.43   46.91 1.00      964     1375
##    m0[24]   45.94   46.06  4.52  4.34   38.24   53.32 1.01      588     1338
##    m0[25]   45.56   45.68  4.43  4.26   38.03   52.77 1.01      605     1303
##    m0[26]   43.54   43.67  3.96  3.73   36.90   49.96 1.01      725     1321
##    m0[27]   20.74   20.78  4.27  4.21   13.37   27.15 1.02      220      243
##    m0[28]   18.02   18.13  4.94  4.91    9.51   25.46 1.02      217      248
##    m0[29]   22.88   22.98  3.78  3.68   16.32   28.57 1.02      227      257
##    m0[30]   25.10   25.18  3.31  3.17   19.26   30.22 1.01      243      299
##    m0[31]   39.11   39.18  3.05  2.83   34.04   44.04 1.00     1262     1357
##    m0[32]   40.25   40.35  3.26  3.10   34.75   45.67 1.00     1092     1404
##    m0[33]   33.81   33.80  2.46  2.38   29.69   37.77 1.00      916     1252
##    m0[34]   22.82   22.91  3.79  3.70   16.25   28.52 1.02      227      257
##    m0[35]   41.80   41.94  3.57  3.37   35.76   47.59 1.00      894     1363
##    m0[36]   23.87   23.96  3.56  3.42   17.58   29.23 1.02      232      275
##    m0[37]   35.74   35.70  2.58  2.48   31.41   39.98 1.00     1437     1347
##    m0[38]   43.69   43.81  3.99  3.78   36.98   50.17 1.01      714     1321
##    m0[39]   28.33   28.35  2.76  2.68   23.58   32.67 1.01      302      458
##    m0[40]   31.64   31.66  2.46  2.41   27.48   35.63 1.00      522      827
##    m0[41]   14.88   14.91  5.75  5.70    4.93   23.57 1.02      216      256
##    m0[42]   42.35   42.50  3.69  3.48   36.12   48.27 1.01      834     1325
##    m0[43]   22.01   22.09  3.97  3.88   15.21   27.98 1.02      224      245
##    m0[44]   25.44   25.51  3.24  3.13   19.78   30.48 1.01      246      301
##    m0[45]   37.80   37.81  2.84  2.64   32.96   42.42 1.00     1443     1326
##    m0[46]   40.67   40.77  3.34  3.19   35.04   46.18 1.00     1037     1402
##    m0[47]   40.90   41.00  3.39  3.23   35.19   46.45 1.00     1008     1401
##    m0[48]   26.52   26.58  3.05  2.91   21.30   31.32 1.01      261      320
##    m0[49]   34.21   34.20  2.47  2.41   30.10   38.20 1.00     1023     1295
##    m0[50]   22.47   22.56  3.87  3.78   15.83   28.30 1.02      225      257
##    y0[1]    19.28   19.53 18.18 17.46  -10.89   48.86 1.00     1528     1827
##    y0[2]    30.58   31.07 17.53 17.53    0.40   58.64 1.00     1891     1957
##    y0[3]    38.38   38.45 17.66 17.84    8.96   66.76 1.00     1851     1913
##    y0[4]    24.63   24.80 17.78 17.50   -5.25   52.96 1.00     1793     2087
##    y0[5]    15.20   15.50 18.09 17.53  -15.41   44.97 1.00     1297     1754
##    y0[6]    18.21   18.47 18.41 17.62  -12.14   47.98 1.00     1573     1357
##    y0[7]    36.38   35.88 17.32 16.86    8.85   65.82 1.00     1681     1895
##    y0[8]    37.63   37.23 17.38 17.22    9.93   66.49 1.00     1879     1869
##    y0[9]    31.34   31.12 17.93 17.53    2.66   60.82 1.00     1833     1933
##    y0[10]   39.58   39.49 18.03 17.82   10.26   69.25 1.00     1977     1833
##    y0[11]   46.69   46.58 18.16 17.87   15.53   75.59 1.00     1481     1368
##    y0[12]   33.40   33.45 17.94 17.88    4.41   62.14 1.00     1721     1973
##    y0[13]   41.98   41.97 17.53 17.23   12.71   71.17 1.00     2064     1970
##    y0[14]   31.77   31.71 17.63 17.42    3.45   60.27 1.00     2199     1985
##    y0[15]   25.20   25.68 17.63 17.21   -4.45   53.96 1.00     1979     1701
##    y0[16]   24.97   25.30 17.57 17.19   -5.30   53.80 1.00     1220     1950
##    y0[17]   43.53   43.43 18.13 17.78   13.37   72.74 1.00     1853     1747
##    y0[18]   32.37   32.25 17.40 16.85    3.99   60.75 1.00     2067     1953
##    y0[19]   41.41   41.05 17.68 17.49   13.66   70.50 1.00     1809     1834
##    y0[20]   35.32   34.97 17.52 16.71    6.83   63.94 1.00     2012     1895
##    y0[21]   29.70   29.61 17.39 17.29    1.20   57.82 1.00     2040     1750
##    y0[22]   20.71   20.86 17.93 17.54   -9.24   49.87 1.00     1532     1971
##    y0[23]   40.82   40.42 17.82 17.01   12.18   70.62 1.00     1796     1562
##    y0[24]   45.77   46.12 18.00 17.66   16.56   76.43 1.00     1761     1836
##    y0[25]   46.14   45.97 18.21 17.97   16.69   76.09 1.00     1846     1930
##    y0[26]   43.59   43.90 17.77 17.39   14.22   72.23 1.00     1723     1871
##    y0[27]   21.05   21.27 18.16 18.64   -8.34   50.24 1.00     1418     1892
##    y0[28]   18.12   18.38 18.66 17.87  -12.75   49.07 1.00     1235     1857
##    y0[29]   21.88   21.83 18.07 17.69   -7.16   51.98 1.00     1488     1966
##    y0[30]   24.81   25.07 17.51 17.30   -5.84   53.77 1.00     1736     1736
##    y0[31]   38.87   38.76 17.40 16.60   10.34   66.48 1.00     1770     1766
##    y0[32]   39.88   39.49 17.59 17.99   11.82   68.76 1.00     2175     2015
##    y0[33]   34.05   34.10 17.65 17.53    4.83   63.03 1.00     2083     1908
##    y0[34]   22.93   22.66 17.93 17.80   -6.49   52.50 1.00     2074     1895
##    y0[35]   41.96   42.07 17.96 17.70   13.19   71.52 1.00     2087     1970
##    y0[36]   23.59   23.31 17.85 17.67   -5.25   53.40 1.00     1818     1998
##    y0[37]   35.62   35.67 17.84 17.83    6.06   64.08 1.00     1961     2012
##    y0[38]   43.46   42.78 17.98 17.73   15.21   74.84 1.00     2023     1674
##    y0[39]   28.20   27.80 17.19 17.86    0.89   57.34 1.00     1863     2004
##    y0[40]   30.84   30.69 17.71 17.80    2.75   59.08 1.00     1978     1930
##    y0[41]   15.19   15.30 18.31 17.65  -15.93   45.09 1.00     1142     1642
##    y0[42]   42.80   43.05 17.87 17.48   13.60   72.10 1.00     1981     1731
##    y0[43]   22.94   22.43 17.76 17.52   -5.66   53.22 1.00     1785     2053
##    y0[44]   25.70   25.73 17.56 16.95   -2.78   55.40 1.00     1719     1971
##    y0[45]   37.84   38.35 18.08 17.44    8.63   67.26 1.00     1821     1787
##    y0[46]   41.34   41.54 17.82 17.38   12.41   70.54 1.00     1994     1820
##    y0[47]   40.67   41.20 17.29 16.97   12.87   69.87 1.00     1943     1709
##    y0[48]   26.72   27.02 17.90 18.13   -2.21   56.44 1.00     1881     1878
##    y0[49]   34.55   34.64 17.71 17.60    5.14   63.59 1.00     2144     2009
##    y0[50]   22.29   22.64 17.46 17.94   -5.89   51.05 1.00     1804     1788
##    m1[1]    14.88   14.91  5.75  5.70    4.93   23.57 1.02      216      256
##    m1[2]    18.36   18.44  4.85  4.78    9.98   25.68 1.02      217      248
##    m1[3]    21.84   21.90  4.01  3.94   14.96   27.88 1.02      223      245
##    m1[4]    25.31   25.38  3.27  3.15   19.58   30.37 1.01      245      297
##    m1[5]    28.79   28.78  2.70  2.63   24.21   33.01 1.01      318      487
##    m1[6]    32.27   32.26  2.44  2.38   28.18   36.26 1.00      605      860
##    m1[7]    35.74   35.70  2.58  2.48   31.41   39.98 1.00     1437     1347
##    m1[8]    39.22   39.32  3.07  2.85   34.09   44.19 1.00     1243     1357
##    m1[9]    42.70   42.82  3.76  3.57   36.35   48.79 1.01      798     1325
##    m1[10]   46.17   46.30  4.58  4.40   38.39   53.70 1.01      578     1338
##    y1[1]    14.72   14.61 18.42 18.70  -15.27   44.10 1.00     1537     1969
##    y1[2]    18.77   19.14 18.18 17.70  -12.25   49.10 1.00     1229     1843
##    y1[3]    22.28   21.90 18.27 17.93   -7.23   52.59 1.00     1724     1843
##    y1[4]    25.33   25.71 17.97 17.62   -5.22   54.43 1.01     1256     1766
##    y1[5]    29.05   28.70 17.41 17.04    0.81   57.79 1.00     1904     1988
##    y1[6]    31.58   31.67 17.69 17.50    3.61   60.51 1.00     1963     1822
##    y1[7]    35.99   36.43 17.54 16.97    7.30   64.84 1.00     2039     1856
##    y1[8]    38.91   38.95 17.95 17.14    9.43   68.88 1.00     2050     2099
##    y1[9]    42.30   42.42 17.99 18.15   12.97   72.76 1.00     1911     1974
##    y1[10]   46.50   46.16 18.14 17.56   16.45   77.16 1.00     1821     1873

estimate as independent class

all b0l,b1l do not have a distribution  
class c=1~k 

ex12-1.stan

yc~N(b0c+b1c*x,s)

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  vector[K] b0;
  vector[K] b1;
  real<lower=0> s;
}
model{
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-1.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -3454.57 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -98.2382     0.0980673      0.370953           1           1      123    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199      -98.2028    0.00186389    0.00452686           1           1      235    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      211      -98.2028   7.21993e-05    0.00177662       0.956       0.956      247    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -98.20
##    b0[1]      5.71
##    b0[2]     -0.30
##    b0[3]     15.80
##    b0[4]      8.88
##    b0[5]     12.25
##    b0[6]     13.78
##    b0[7]     16.90
##    b0[8]      9.46
##    b0[9]      8.57
##    b1[1]      0.84
##    b1[2]      4.24
##    b1[3]     -0.38
##    b1[4]      4.18
##    b1[5]      0.18
##    b1[6]      2.01
##    b1[7]      2.02
##    b1[8]      2.61
##    b1[9]      2.73
##    s          4.32
##    m0[1]     22.32
##    m0[2]     35.32
##    m0[3]     48.95
##    m0[4]     29.75
##    m0[5]     15.66
##    m0[6]     15.25
##    m0[7]     56.51
##    m0[8]     10.16
##    m0[9]     14.15
##    m0[10]    48.87
##    m0[11]     8.22
##    m0[12]    14.37
##    m0[13]    54.36
##    m0[14]    37.45
##    m0[15]    25.93
##    m0[16]    26.58
##    m0[17]    20.65
##    m0[18]    14.26
##    m0[19]    51.53
##    m0[20]    42.89
##    m0[21]    50.09
##    m0[22]    21.82
##    m0[23]    54.06
##    m0[24]    53.20
##    m0[25]    61.34
##    m0[26]    15.56
##    m0[27]    19.14
##    m0[28]    21.31
##    m0[29]    13.80
##    m0[30]    26.67
##    m0[31]    50.37
##    m0[32]    46.06
##    m0[33]    41.22
##    m0[34]    24.19
##    m0[35]    71.75
##    m0[36]    13.31
##    m0[37]    10.72
##    m0[38]    84.94
##    m0[39]    36.16
##    m0[40]    53.51
##    m0[41]     5.90
##    m0[42]     9.14
##    m0[43]    21.29
##    m0[44]    28.54
##    m0[45]    17.83
##    m0[46]    19.32
##    m0[47]    50.17
##    m0[48]    32.03
##    m0[49]    42.04
##    m0[50]    26.93
##    y0[1]     26.33
##    y0[2]     36.61
##    y0[3]     45.89
##    y0[4]     28.22
##    y0[5]     19.82
##    y0[6]     14.93
##    y0[7]     54.59
##    y0[8]     11.15
##    y0[9]     14.70
##    y0[10]    53.16
##    y0[11]     3.03
##    y0[12]    19.77
##    y0[13]    60.89
##    y0[14]    41.39
##    y0[15]    29.59
##    y0[16]    23.59
##    y0[17]    24.79
##    y0[18]    13.43
##    y0[19]    53.13
##    y0[20]    48.38
##    y0[21]    45.98
##    y0[22]    16.53
##    y0[23]    47.89
##    y0[24]    52.90
##    y0[25]    68.99
##    y0[26]    13.91
##    y0[27]    17.29
##    y0[28]    18.32
##    y0[29]     7.09
##    y0[30]    24.85
##    y0[31]    50.87
##    y0[32]    44.12
##    y0[33]    53.86
##    y0[34]    23.39
##    y0[35]    66.00
##    y0[36]    13.44
##    y0[37]    13.45
##    y0[38]    84.43
##    y0[39]    38.67
##    y0[40]    60.75
##    y0[41]     9.40
##    y0[42]     3.88
##    y0[43]    19.34
##    y0[44]    29.15
##    y0[45]    18.83
##    y0[46]    12.34
##    y0[47]    47.24
##    y0[48]    29.06
##    y0[49]    46.98
##    y0[50]    39.90
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 1.6 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 1.6 seconds.
## Chain 3 finished in 1.6 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 1.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.6 seconds.
## Total execution time: 1.8 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -108.74 -108.29  3.82  3.69 -115.91 -103.30 1.00      570      818
##    b0[1]     5.79    5.77  5.54  5.38   -3.22   14.84 1.01     2301      937
##    b0[2]    -0.36   -0.64  8.67  8.56  -14.09   14.14 1.00     1270     1262
##    b0[3]    15.79   15.83  4.64  4.51    8.35   23.54 1.00     3237     1460
##    b0[4]     9.04    8.55 11.58 11.35   -8.82   28.65 1.01      713      829
##    b0[5]    12.46   12.52  7.68  7.97   -0.30   24.84 1.00     1328     1402
##    b0[6]    13.84   13.90  5.56  5.37    4.75   23.10 1.00     3017     1255
##    b0[7]    16.94   16.89  3.72  3.53   10.71   23.02 1.00     3335     1377
##    b0[8]     9.79    9.65  7.31  7.14   -2.12   22.15 1.00     1840     1355
##    b0[9]     8.57    8.53  3.94  3.90    2.15   15.33 1.00     2956      987
##    b1[1]     0.83    0.84  0.39  0.39    0.17    1.45 1.00     2074     1446
##    b1[2]     4.25    4.26  0.71  0.70    3.09    5.39 1.00     1335     1421
##    b1[3]    -0.39   -0.38  0.34  0.34   -0.95    0.17 1.00     2334     1594
##    b1[4]     4.16    4.22  0.86  0.84    2.71    5.48 1.01      723      918
##    b1[5]     0.18    0.17  0.63  0.65   -0.82    1.22 1.00     1360     1553
##    b1[6]     2.01    2.00  0.42  0.40    1.32    2.71 1.00     2784     1626
##    b1[7]     2.03    2.03  0.34  0.33    1.48    2.59 1.00     2861     1417
##    b1[8]     2.59    2.59  0.66  0.65    1.49    3.68 1.00     1904     1333
##    b1[9]     2.73    2.74  0.32  0.32    2.18    3.25 1.00     2676     1590
##    s         5.61    5.57  0.71  0.68    4.57    6.88 1.00     1035     1337
##    m0[1]    22.36   22.38  2.96  2.76   17.47   27.26 1.00     3401     1376
##    m0[2]    35.42   35.37  2.55  2.49   31.36   39.61 1.00     1960     1740
##    m0[3]    48.99   49.07  2.31  2.25   45.11   52.77 1.00     1915     1417
##    m0[4]    29.80   29.81  2.13  2.02   26.23   33.27 1.00     3020     1521
##    m0[5]    15.65   15.67  4.53  4.38    8.37   23.25 1.00     3249     1460
##    m0[6]    15.26   15.18  3.27  3.20    9.79   20.75 1.00     2952     1165
##    m0[7]    56.61   56.56  3.26  3.21   51.40   62.10 1.00     2003     1656
##    m0[8]    10.12   10.07  2.58  2.62    5.98   14.33 1.00     1618     1149
##    m0[9]    14.29   14.33  2.54  2.54   10.08   18.43 1.00     1740     1531
##    m0[10]   48.95   48.96  2.77  2.73   44.33   53.44 1.01     2130     1365
##    m0[11]    8.16    8.20  3.65  3.66    2.19   14.17 1.00     1585     1619
##    m0[12]   14.50   14.54  2.44  2.43   10.49   18.51 1.00     1863     1505
##    m0[13]   54.30   54.17  5.17  5.08   46.07   62.61 1.00     2046     1500
##    m0[14]   37.54   37.55  2.54  2.47   33.51   41.74 1.00     1953     1698
##    m0[15]   25.95   25.98  2.36  2.33   22.10   29.76 1.00     2841     1377
##    m0[16]   26.75   26.72  3.57  3.49   20.89   32.87 1.00     1838     1443
##    m0[17]   20.63   20.66  3.60  3.40   14.75   26.61 1.00     1967     1418
##    m0[18]   14.39   14.42  2.46  2.47   10.31   18.47 1.00     1804     1439
##    m0[19]   51.61   51.56  3.12  3.13   46.49   56.70 1.01     2186     1429
##    m0[20]   42.96   42.98  2.13  2.08   39.42   46.40 1.01     2056     1167
##    m0[21]   50.05   50.00  4.19  4.03   43.40   57.16 1.00      946     1146
##    m0[22]   21.87   21.77  4.20  4.08   15.10   28.91 1.00     2983     1140
##    m0[23]   54.11   54.16  2.70  2.54   49.54   58.56 1.01     1989     1338
##    m0[24]   53.17   53.21  4.38  4.33   45.96   60.45 1.00     2113     1677
##    m0[25]   61.39   61.43  3.36  3.22   55.76   66.93 1.01     2124     1551
##    m0[26]   15.65   15.64  4.83  4.75    7.88   23.49 1.00     1643     1540
##    m0[27]   19.15   19.09  2.91  2.87   14.30   23.95 1.00     2924     1260
##    m0[28]   21.35   21.37  3.10  2.87   16.21   26.44 1.00     3423     1377
##    m0[29]   13.78   13.81  3.22  3.16    8.60   19.03 1.00     3110     1356
##    m0[30]   26.84   26.80  3.56  3.48   20.98   32.93 1.00     1839     1421
##    m0[31]   50.41   50.47  2.41  2.30   46.32   54.33 1.00     1931     1379
##    m0[32]   46.05   46.07  3.37  3.29   40.47   51.64 1.00     2046     1445
##    m0[33]   41.29   41.33  2.00  1.96   38.03   44.53 1.01     2081     1294
##    m0[34]   24.23   24.14  3.85  3.77   18.01   30.61 1.00     2943     1138
##    m0[35]   71.90   71.99  4.94  4.98   63.60   80.01 1.00     1731     1487
##    m0[36]   13.49   13.54  4.34  4.56    6.25   20.53 1.00     1356     1575
##    m0[37]   10.67   10.68  2.42  2.47    6.80   14.64 1.00     1701     1366
##    m0[38]   84.73   84.90  5.51  5.21   75.47   93.44 1.00     1106     1345
##    m0[39]   36.20   36.21  3.56  3.48   30.27   42.20 1.00     1464     1422
##    m0[40]   53.45   53.37  3.78  3.66   47.40   59.87 1.00     1075     1224
##    m0[41]    5.98    5.97  5.46  5.36   -2.92   14.92 1.01     2304      937
##    m0[42]    9.08    9.12  3.07  3.11    4.15   14.12 1.00     1576     1291
##    m0[43]   21.31   21.27  2.73  2.70   16.78   25.79 1.00     2893     1259
##    m0[44]   28.56   28.51  4.42  4.28   21.40   36.00 1.00     1341     1331
##    m0[45]   17.84   17.83  2.94  2.85   12.98   22.74 1.00     2003     1554
##    m0[46]   19.32   19.34  3.24  3.05   13.95   24.70 1.00     1970     1573
##    m0[47]   50.24   50.21  2.94  2.91   45.40   55.02 1.01     2160     1338
##    m0[48]   32.08   32.08  1.96  1.85   28.88   35.18 1.00     2826     1432
##    m0[49]   42.07   42.05  1.96  1.86   38.92   45.20 1.00     1975     1293
##    m0[50]   26.97   27.00  2.41  2.30   22.95   30.93 1.00     3145     1520
##    y0[1]    22.38   22.34  6.18  6.04   12.18   32.72 1.00     2115     2055
##    y0[2]    35.52   35.66  5.98  5.99   25.74   45.32 1.00     2192     1972
##    y0[3]    48.90   49.10  6.06  5.99   38.75   58.79 1.00     2036     1996
##    y0[4]    29.67   29.86  6.19  6.06   19.22   40.16 1.00     2252     1845
##    y0[5]    15.46   15.48  7.31  7.09    3.06   27.56 1.00     2363     1685
##    y0[6]    15.18   15.26  6.67  6.47    4.19   26.03 1.00     2347     2098
##    y0[7]    56.52   56.51  6.52  6.53   45.52   67.42 1.00     1943     1892
##    y0[8]    10.01   10.25  6.25  6.15   -0.11   20.33 1.00     1968     1916
##    y0[9]    14.16   14.19  6.20  6.20    4.18   24.50 1.00     2091     1952
##    y0[10]   49.08   49.06  6.46  6.45   38.57   59.71 1.00     2116     1598
##    y0[11]    8.31    8.21  6.75  6.76   -2.58   19.51 1.00     1930     1787
##    y0[12]   14.46   14.45  6.08  5.90    4.31   24.48 1.00     1898     1881
##    y0[13]   54.28   54.35  7.82  7.59   41.47   66.65 1.00     1883     1772
##    y0[14]   37.35   37.48  6.26  5.88   26.89   47.68 1.00     2057     2054
##    y0[15]   25.70   25.62  6.02  5.79   15.92   35.56 1.00     1892     1726
##    y0[16]   26.68   26.71  6.78  6.81   15.42   37.30 1.00     1877     1836
##    y0[17]   20.49   20.31  6.88  6.88    9.46   31.83 1.00     2228     1748
##    y0[18]   14.44   14.40  6.22  6.11    4.08   24.64 1.00     1902     1688
##    y0[19]   51.50   51.69  6.47  6.41   40.97   62.02 1.00     1776     1634
##    y0[20]   43.05   42.99  6.00  5.67   33.24   52.97 1.00     2035     1892
##    y0[21]   50.21   50.29  7.02  6.99   38.80   61.73 1.00     1313     1582
##    y0[22]   21.77   21.89  7.12  7.06    9.92   33.11 1.00     2380     1892
##    y0[23]   54.15   54.35  6.24  6.21   43.32   63.75 1.00     1928     1788
##    y0[24]   53.16   53.17  7.18  7.09   41.58   65.32 1.00     2069     1930
##    y0[25]   61.57   61.58  6.54  6.30   50.90   72.28 1.00     2110     1894
##    y0[26]   15.63   15.63  7.30  7.40    3.91   28.20 1.00     1900     1857
##    y0[27]   19.06   19.01  6.25  5.94    8.88   29.03 1.00     2171     1999
##    y0[28]   21.65   21.47  6.56  6.41   10.98   32.60 1.00     2220     1971
##    y0[29]   13.76   13.54  6.51  6.38    3.53   24.65 1.00     1853     1791
##    y0[30]   26.77   26.70  6.72  6.43   15.55   37.89 1.00     1978     1829
##    y0[31]   50.45   50.46  6.25  6.06   40.09   61.17 1.00     2020     1923
##    y0[32]   46.07   46.13  6.56  6.54   35.24   56.67 1.00     2031     1832
##    y0[33]   41.26   41.35  6.08  5.95   31.36   51.02 1.00     2128     1931
##    y0[34]   23.93   23.86  6.85  6.65   12.88   35.06 1.00     2397     1905
##    y0[35]   71.93   71.92  7.64  7.50   59.56   84.53 1.00     1887     1880
##    y0[36]   13.43   13.53  7.03  7.00    2.18   24.79 1.00     1770     1852
##    y0[37]   10.59   10.47  6.12  6.15    0.60   20.92 1.00     1976     2014
##    y0[38]   84.78   85.15  7.91  7.73   71.60   97.20 1.00     1487     1668
##    y0[39]   36.30   36.25  6.72  6.44   25.36   47.29 1.00     1985     1931
##    y0[40]   53.30   53.23  6.85  6.56   41.89   64.80 1.00     1512     1708
##    y0[41]    5.99    6.15  7.93  7.94   -6.78   18.88 1.00     2044     1710
##    y0[42]    9.25    9.28  6.38  6.26   -1.22   19.44 1.00     1570     1879
##    y0[43]   21.42   21.22  6.18  6.17   11.81   31.56 1.00     2062     1974
##    y0[44]   28.47   28.41  7.28  7.19   16.51   41.03 1.00     1809     1819
##    y0[45]   17.92   17.92  6.30  6.18    7.54   28.40 1.00     1968     1799
##    y0[46]   19.24   19.12  6.46  6.35    8.51   29.99 1.00     1904     1750
##    y0[47]   50.53   50.43  6.31  6.29   40.28   61.10 1.00     2129     1449
##    y0[48]   31.88   31.88  5.85  5.83   21.97   41.14 1.00     1936     1777
##    y0[49]   41.97   41.91  5.89  5.79   32.24   51.91 1.00     2096     2057
##    y0[50]   26.79   26.74  6.02  6.07   17.09   36.57 1.00     2015     2014

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

estimate as class have relation

all b0l,b1l have a distribution  
class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)

ex12-2.stan

class have relation

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  real b00;
  real<lower=0,upper=100> sb0;
  vector[K] b0;
  real b10;
  real<lower=0,upper=100> sb1;
  vector[K] b1;
  real<lower=0,upper=100> s;
}
model{
  b0~normal(b00,sb0);
  b1~normal(b10,sb1);
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-2.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -298.118 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -29.1277      0.280031   1.12786e+06      0.6334      0.6334      166    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      163     0.0962256   3.18287e-05   1.38466e+06      0.8781      0.8781      256    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
try(print(mle))
##  variable estimate
##    lp__       0.10
##    b00        0.92
##    sb0        0.00
##    b0[1]      0.92
##    b0[2]      0.92
##    b0[3]      0.92
##    b0[4]      0.92
##    b0[5]      0.92
##    b0[6]      0.92
##    b0[7]      0.92
##    b0[8]      0.92
##    b0[9]      0.92
##    b10        2.14
##    sb1        1.07
##    b1[1]      1.88
##    b1[2]      4.59
##    b1[3]      1.92
##    b1[4]      4.38
##    b1[5]      2.30
##    b1[6]      1.67
##    b1[7]      2.94
##    b1[8]      2.73
##    b1[9]      2.23
##    s          6.31
##    m0[1]      8.81
##    m0[2]     27.98
##    m0[3]     33.94
##    m0[4]     19.61
##    m0[5]      1.62
##    m0[6]      6.39
##    m0[7]     62.47
##    m0[8]     29.01
##    m0[9]     24.78
##    m0[10]    47.43
##    m0[11]    38.71
##    m0[12]    27.55
##    m0[13]    47.92
##    m0[14]    30.22
##    m0[15]    15.12
##    m0[16]    18.83
##    m0[17]    34.60
##    m0[18]    26.14
##    m0[19]    51.29
##    m0[20]    38.73
##    m0[21]    44.10
##    m0[22]     7.61
##    m0[23]    38.12
##    m0[24]    33.72
##    m0[25]    44.06
##    m0[26]    42.55
##    m0[27]     9.56
##    m0[28]     7.33
##    m0[29]    10.90
##    m0[30]    18.93
##    m0[31]    35.10
##    m0[32]    27.78
##    m0[33]    36.30
##    m0[34]     9.58
##    m0[35]    78.98
##    m0[36]    14.33
##    m0[37]    26.26
##    m0[38]    80.61
##    m0[39]    40.42
##    m0[40]    47.68
##    m0[41]     1.35
##    m0[42]    34.14
##    m0[43]    11.33
##    m0[44]    32.17
##    m0[45]    28.25
##    m0[46]    31.62
##    m0[47]    49.31
##    m0[48]    22.93
##    m0[49]    28.28
##    m0[50]    15.51
##    y0[1]     13.94
##    y0[2]     30.95
##    y0[3]     36.00
##    y0[4]     18.97
##    y0[5]      6.09
##    y0[6]     10.41
##    y0[7]     63.29
##    y0[8]     23.04
##    y0[9]     21.54
##    y0[10]    61.53
##    y0[11]    36.39
##    y0[12]    21.56
##    y0[13]    46.39
##    y0[14]    37.40
##    y0[15]    16.80
##    y0[16]    16.12
##    y0[17]    24.91
##    y0[18]    26.44
##    y0[19]    56.45
##    y0[20]    37.97
##    y0[21]    46.34
##    y0[22]     4.63
##    y0[23]    38.75
##    y0[24]    43.31
##    y0[25]    44.43
##    y0[26]    45.71
##    y0[27]    15.61
##    y0[28]    12.90
##    y0[29]     9.94
##    y0[30]    14.84
##    y0[31]    22.51
##    y0[32]    33.20
##    y0[33]    37.90
##    y0[34]    22.00
##    y0[35]    80.42
##    y0[36]     7.74
##    y0[37]    19.19
##    y0[38]    84.07
##    y0[39]    33.08
##    y0[40]    36.25
##    y0[41]    10.50
##    y0[42]    40.71
##    y0[43]    14.35
##    y0[44]    27.35
##    y0[45]    27.76
##    y0[46]    21.52
##    y0[47]    38.06
##    y0[48]    20.08
##    y0[49]    27.12
##    y0[50]    27.27
mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 3 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 3 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 4 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 4 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 3 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 4 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 1 finished in 1.3 seconds.
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 3 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 4 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 finished in 1.4 seconds.
## Chain 3 finished in 1.3 seconds.
## Chain 4 finished in 1.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.4 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -127.67 -127.82 6.45 6.14 -138.27 -116.43 1.00      484      284
##    b00      11.51   11.57 2.36 2.20    7.67   15.25 1.00     2577     3014
##    sb0       3.58    3.04 2.47 2.13    0.75    8.24 1.00      525      273
##    b0[1]     9.66   10.03 3.62 3.30    3.04   14.89 1.00     2747     4324
##    b0[2]    10.11   10.48 4.07 3.47    2.75   15.98 1.00     3137     4413
##    b0[3]    12.64   12.52 3.03 2.80    7.87   17.75 1.00     4199     4841
##    b0[4]    12.06   11.97 4.24 3.39    5.59   19.13 1.00     4565     3786
##    b0[5]    11.19   11.34 3.74 3.32    4.86   16.87 1.00     3977     4723
##    b0[6]    12.29   12.19 3.19 2.94    7.30   17.56 1.00     4372     4316
##    b0[7]    13.99   13.75 2.97 2.82    9.45   19.19 1.00     2132     2293
##    b0[8]    11.29   11.44 3.49 3.06    5.32   16.82 1.00     4323     4076
##    b0[9]    10.51   10.64 2.84 2.69    5.60   14.92 1.00     3084     3632
##    b10       1.95    1.95 0.61 0.54    0.97    2.92 1.00     4096     2949
##    sb1       1.70    1.59 0.57 0.45    1.04    2.74 1.00     2302     1311
##    b1[1]     0.61    0.60 0.29 0.28    0.13    1.10 1.00     2699      940
##    b1[2]     3.39    3.38 0.39 0.36    2.78    4.06 1.00     3045     1514
##    b1[3]    -0.16   -0.15 0.26 0.25   -0.58    0.27 1.00     2179     1097
##    b1[4]     3.90    3.90 0.38 0.35    3.29    4.52 1.00     3926     3895
##    b1[5]     0.30    0.29 0.36 0.33   -0.27    0.88 1.00     4740     4861
##    b1[6]     2.11    2.11 0.29 0.28    1.63    2.58 1.00     4518     5289
##    b1[7]     2.25    2.26 0.28 0.26    1.77    2.68 1.00     2342     2776
##    b1[8]     2.44    2.43 0.37 0.34    1.83    3.04 1.00     4995     4857
##    b1[9]     2.58    2.58 0.25 0.25    2.18    3.00 1.00     3815     4565
##    s         5.47    5.40 0.67 0.64    4.50    6.70 1.00     2562     1449
##    m0[1]    20.01   19.85 2.39 2.28   16.35   24.19 1.00     2343     3032
##    m0[2]    35.42   35.44 2.25 2.18   31.71   39.15 1.00     5823     4865
##    m0[3]    48.74   48.72 2.21 2.15   45.07   52.39 1.00     6963     4519
##    m0[4]    28.26   28.19 1.80 1.73   25.44   31.32 1.00     3449     5194
##    m0[5]    12.58   12.46 2.96 2.74    7.87   17.59 1.00     4227     4791
##    m0[6]    16.84   16.95 2.38 2.23   12.75   20.48 1.00     3135     3448
##    m0[7]    55.57   55.58 3.04 3.04   50.54   60.57 1.00     5032     1637
##    m0[8]    10.32   10.33 2.54 2.46    6.18   14.44 1.00     2105      966
##    m0[9]    14.27   14.26 2.22 2.15   10.63   17.89 1.00     5654     4662
##    m0[10]   49.50   49.52 2.53 2.50   45.37   53.60 1.00     5780     5241
##    m0[11]    9.51    9.51 3.49 3.34    3.89   15.25 1.00     1879      952
##    m0[12]   14.62   14.61 2.37 2.27   10.75   18.52 1.00     6008     4650
##    m0[13]   53.20   53.23 4.11 4.02   46.54   60.01 1.00     5880     5218
##    m0[14]   37.42   37.44 2.38 2.28   33.50   41.32 1.00     5962     5026
##    m0[15]   26.95   27.01 1.82 1.77   23.87   29.80 1.00     3694     3905
##    m0[16]   27.26   27.29 2.12 2.04   23.78   30.65 1.00     4930     4995
##    m0[17]   20.55   20.50 3.48 3.46   14.94   26.25 1.00     3457     3332
##    m0[18]   14.44   14.43 2.28 2.21   10.69   18.19 1.00     5866     4765
##    m0[19]   52.45   52.49 2.82 2.76   47.82   57.01 1.00     5016     5446
##    m0[20]   42.86   42.86 1.99 1.90   39.60   46.08 1.00     8074     5164
##    m0[21]   50.50   50.55 2.68 2.62   46.01   54.84 1.00     2976     1030
##    m0[22]   20.72   20.64 2.51 2.27   16.76   24.95 1.00     4860     4403
##    m0[23]   53.59   53.55 2.54 2.50   49.38   57.76 1.00     6697     4719
##    m0[24]   53.63   53.65 4.18 4.04   46.73   60.47 1.00     5691     5324
##    m0[25]   60.47   60.47 3.08 2.99   55.38   65.47 1.00     6116     6405
##    m0[26]   16.56   16.46 3.98 3.88   10.06   23.18 1.00     5725     5864
##    m0[27]   20.51   20.61 2.14 2.03   16.81   23.77 1.00     3235     3792
##    m0[28]   18.89   18.70 2.50 2.37   15.06   23.22 1.00     2285     3081
##    m0[29]   11.82   11.73 2.22 2.11    8.27   15.59 1.00     4459     4878
##    m0[30]   27.35   27.38 2.12 2.05   23.86   30.74 1.00     4937     5011
##    m0[31]   50.09   50.06 2.30 2.24   46.28   53.89 1.00     6901     4690
##    m0[32]   46.15   46.12 3.36 3.21   40.54   51.68 1.00     5979     5008
##    m0[33]   41.00   41.00 1.87 1.78   37.93   44.04 1.00     8174     5268
##    m0[34]   23.20   23.13 2.39 2.20   19.41   27.24 1.00     5055     4138
##    m0[35]   67.77   67.85 4.08 4.01   60.90   74.43 1.00     4162     1627
##    m0[36]   12.92   12.96 2.36 2.17    8.93   16.60 1.00     4153     5037
##    m0[37]   10.54   10.56 2.34 2.26    6.75   14.36 1.00     2252      967
##    m0[38]   83.00   83.10 4.51 4.38   75.20   90.35 1.00     2951      976
##    m0[39]   39.29   39.30 2.36 2.30   35.43   43.09 1.00     5058     5867
##    m0[40]   53.68   53.75 2.75 2.69   49.07   58.16 1.00     2798     1018
##    m0[41]    9.80   10.16 3.57 3.27    3.26   14.96 1.00     2758     4443
##    m0[42]    9.89    9.86 3.01 2.92    5.00   14.76 1.00     1933      988
##    m0[43]   22.56   22.65 2.02 1.93   19.08   25.67 1.00     3335     3958
##    m0[44]   33.19   33.25 2.45 2.30   29.05   36.98 1.00     4399     5433
##    m0[45]   18.50   18.45 2.86 2.81   13.90   23.30 1.00     3973     4705
##    m0[46]   19.59   19.53 3.17 3.12   14.48   24.89 1.00     3689     4840
##    m0[47]   50.94   50.97 2.67 2.63   46.55   55.26 1.00     5367     5460
##    m0[48]   30.79   30.74 1.71 1.64   28.09   33.68 1.00     4333     5527
##    m0[49]   42.20   42.20 1.86 1.81   39.13   45.27 1.00     6548     4363
##    m0[50]   25.13   25.02 1.99 1.92   22.05   28.52 1.00     2840     4361
##    y0[1]    19.94   19.87 6.05 5.85   10.00   29.94 1.00     6955     7166
##    y0[2]    35.36   35.37 5.91 5.78   25.65   45.00 1.00     7888     7091
##    y0[3]    48.69   48.74 6.00 5.95   38.90   58.53 1.00     8348     7735
##    y0[4]    28.16   28.20 5.76 5.73   18.66   37.58 1.00     7904     7642
##    y0[5]    12.65   12.70 6.34 6.26    2.21   22.90 1.00     6808     7459
##    y0[6]    16.81   16.88 6.02 5.90    7.04   26.63 1.00     7310     7592
##    y0[7]    55.52   55.51 6.24 6.04   45.16   65.61 1.00     8259     7526
##    y0[8]    10.31   10.35 6.07 5.96    0.30   20.27 1.00     5976     7655
##    y0[9]    14.33   14.36 5.88 5.79    4.61   23.92 1.00     8158     7779
##    y0[10]   49.55   49.61 6.05 5.95   39.67   59.39 1.00     7730     7677
##    y0[11]    9.44    9.46 6.49 6.42   -1.01   19.91 1.00     4661     4641
##    y0[12]   14.57   14.53 5.98 5.84    4.94   24.52 1.00     7647     7719
##    y0[13]   53.24   53.32 6.93 6.68   41.78   64.41 1.00     7089     7608
##    y0[14]   37.44   37.45 5.96 5.86   27.78   47.22 1.00     7957     7632
##    y0[15]   27.00   26.95 5.84 5.85   17.61   36.47 1.00     7640     7800
##    y0[16]   27.27   27.25 5.82 5.80   17.68   36.79 1.00     7342     7930
##    y0[17]   20.62   20.62 6.52 6.25    9.84   31.38 1.00     7453     7719
##    y0[18]   14.36   14.45 6.01 5.84    4.24   24.25 1.00     8288     8057
##    y0[19]   52.47   52.39 6.13 5.96   42.68   62.60 1.00     7836     7847
##    y0[20]   42.88   42.86 5.84 5.75   33.25   52.41 1.00     8276     7855
##    y0[21]   50.44   50.39 6.14 5.98   40.32   60.51 1.00     7520     7321
##    y0[22]   20.79   20.72 6.01 5.92   11.10   30.75 1.00     7117     7280
##    y0[23]   53.57   53.56 6.00 5.93   43.74   63.24 1.00     7770     6694
##    y0[24]   53.62   53.52 6.87 6.70   42.32   64.88 1.00     7546     6710
##    y0[25]   60.51   60.48 6.30 6.11   50.26   70.76 1.00     8320     7630
##    y0[26]   16.59   16.56 6.81 6.63    5.52   28.03 1.00     7209     7013
##    y0[27]   20.42   20.46 5.88 5.77   10.69   30.12 1.00     6954     7102
##    y0[28]   18.96   18.99 6.05 6.08    9.06   28.82 1.00     6721     7176
##    y0[29]   11.80   11.83 5.92 5.94    2.21   21.49 1.00     8079     7972
##    y0[30]   27.27   27.28 5.89 5.73   17.52   36.94 1.00     7868     7731
##    y0[31]   50.16   50.14 5.96 5.82   40.48   60.14 1.00     7698     7812
##    y0[32]   46.14   46.22 6.41 6.37   35.83   56.73 1.00     7767     7435
##    y0[33]   40.94   40.98 5.80 5.63   31.35   50.45 1.00     7887     7649
##    y0[34]   23.12   23.10 6.06 5.93   13.34   33.13 1.00     7822     7161
##    y0[35]   67.74   67.83 6.82 6.58   56.45   78.81 1.00     5690     6428
##    y0[36]   13.00   13.01 5.98 5.84    3.10   22.81 1.00     6983     7602
##    y0[37]   10.60   10.65 6.02 5.84    0.66   20.48 1.00     6856     7218
##    y0[38]   83.02   83.07 7.16 7.07   71.14   94.68 1.00     4983     5998
##    y0[39]   39.34   39.39 5.93 5.94   29.79   49.02 1.00     7688     7518
##    y0[40]   53.63   53.68 6.19 6.14   43.34   63.67 1.00     6851     7788
##    y0[41]    9.81    9.72 6.57 6.50   -0.93   20.56 1.00     4819     6683
##    y0[42]    9.83    9.81 6.30 6.13   -0.45   20.23 1.00     5210     6413
##    y0[43]   22.58   22.57 5.83 5.75   13.10   32.15 1.00     6605     7333
##    y0[44]   33.30   33.31 6.09 6.03   23.29   43.33 1.00     7359     7834
##    y0[45]   18.43   18.54 6.25 6.19    7.97   28.50 1.00     7375     6806
##    y0[46]   19.56   19.59 6.34 6.26    9.13   30.03 1.00     6599     7050
##    y0[47]   50.95   50.97 6.15 6.14   41.09   61.11 1.00     7560     7443
##    y0[48]   30.83   30.88 5.74 5.75   21.33   40.17 1.00     7571     7651
##    y0[49]   42.18   42.22 5.76 5.67   32.74   51.87 1.00     7699     7161
##    y0[50]   25.16   25.14 5.86 5.84   15.49   34.87 1.00     6799     7173

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

generalized linear mixed model

(X,y)i=1-n
b[b0,b1,...]
linear model    y~N(Xb,s)  
generalized linear model    y~dist.(m=link(Xb),s)  

fixed effect    b0, b1,...  
individual random effect   b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...  
class c=1-k  
class effect    b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...  

for y=dist.(m,s)
random intercept model    m=b0+r0i+b1*x, m=b0+r0c+b1*x  
random coefficient model  m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x  
mixed model   m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x  

note  
@ yi=b0+b1*xi+r0i is not useful, ri is included in s  
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
  but variance is larger than mean (over dispersion),
  random intercept model is useful
n=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)

ex13-0.stan

generalized linear model, poisson regression

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~poisson_log(b0+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-0.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -3.4008e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       18       495.512    0.00132099     0.0675192           1           1       48    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     495.51
##    b0         2.23
##    b1         0.09
##    m0[1]      2.24
##    m0[2]      2.25
##    m0[3]      2.44
##    m0[4]      2.30
##    m0[5]      3.01
##    m0[6]      2.82
##    m0[7]      2.52
##    m0[8]      2.94
##    m0[9]      3.05
##    m0[10]     2.76
##    m0[11]     2.62
##    m0[12]     2.52
##    m0[13]     3.02
##    m0[14]     2.95
##    m0[15]     2.66
##    m0[16]     2.55
##    m0[17]     2.38
##    m0[18]     2.54
##    m0[19]     2.86
##    m0[20]     2.35
##    y0[1]      9.00
##    y0[2]      8.00
##    y0[3]     11.00
##    y0[4]     14.00
##    y0[5]     19.00
##    y0[6]     17.00
##    y0[7]      6.00
##    y0[8]     11.00
##    y0[9]     20.00
##    y0[10]    18.00
##    y0[11]    17.00
##    y0[12]    13.00
##    y0[13]    18.00
##    y0[14]    20.00
##    y0[15]    12.00
##    y0[16]    18.00
##    y0[17]    10.00
##    y0[18]    17.00
##    y0[19]    13.00
##    y0[20]    14.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   494.48 494.83 1.02 0.73 492.41 495.46 1.00      639      984
##    b0       2.23   2.23 0.12 0.12   2.02   2.43 1.01      380      404
##    b1       0.09   0.09 0.02 0.02   0.05   0.12 1.00      435      655
##    m0[1]    2.24   2.24 0.12 0.12   2.03   2.44 1.01      381      404
##    m0[2]    2.25   2.25 0.12 0.12   2.05   2.45 1.01      381      404
##    m0[3]    2.44   2.44 0.08 0.08   2.30   2.58 1.01      416      495
##    m0[4]    2.30   2.30 0.11 0.11   2.12   2.48 1.01      384      417
##    m0[5]    3.00   3.00 0.09 0.09   2.85   3.16 1.00     1129     1200
##    m0[6]    2.82   2.82 0.07 0.06   2.72   2.93 1.00     1767     1477
##    m0[7]    2.51   2.51 0.07 0.07   2.39   2.63 1.01      463      546
##    m0[8]    2.94   2.94 0.08 0.08   2.80   3.08 1.00     1426     1332
##    m0[9]    3.05   3.05 0.10 0.10   2.88   3.22 1.00     1009     1070
##    m0[10]   2.76   2.76 0.06 0.06   2.66   2.86 1.00     1536     1302
##    m0[11]   2.62   2.62 0.06 0.06   2.51   2.72 1.00      654      897
##    m0[12]   2.51   2.51 0.07 0.07   2.39   2.63 1.01      463      546
##    m0[13]   3.01   3.01 0.09 0.09   2.86   3.17 1.00     1104     1200
##    m0[14]   2.95   2.95 0.08 0.08   2.81   3.09 1.00     1371     1351
##    m0[15]   2.66   2.65 0.06 0.06   2.55   2.76 1.00      809      987
##    m0[16]   2.55   2.55 0.07 0.07   2.44   2.67 1.00      511      720
##    m0[17]   2.38   2.38 0.09 0.09   2.22   2.54 1.01      396      433
##    m0[18]   2.54   2.54 0.07 0.07   2.42   2.66 1.00      493      650
##    m0[19]   2.85   2.85 0.07 0.07   2.74   2.97 1.00     1742     1441
##    m0[20]   2.35   2.35 0.10 0.10   2.18   2.51 1.01      390      423
##    y0[1]    9.54   9.00 3.32 2.97   5.00  15.00 1.00     1279     1882
##    y0[2]    9.56   9.00 3.31 2.97   5.00  16.00 1.00     1506     1568
##    y0[3]   11.47  11.00 3.52 2.97   6.00  18.00 1.00     1478     1928
##    y0[4]    9.99  10.00 3.29 2.97   5.00  16.00 1.00     1655     1773
##    y0[5]   20.28  20.00 4.86 4.45  13.00  29.00 1.00     1644     1847
##    y0[6]   16.72  17.00 4.23 4.45  10.00  24.00 1.00     1738     1880
##    y0[7]   12.47  12.00 3.74 4.45   7.00  19.00 1.00     1741     1761
##    y0[8]   18.94  19.00 4.66 4.45  12.00  27.00 1.00     1852     1945
##    y0[9]   21.03  21.00 5.01 4.45  13.00  30.00 1.00     1803     1839
##    y0[10]  15.76  16.00 4.12 4.45   9.00  23.00 1.00     1776     1845
##    y0[11]  13.85  14.00 3.77 2.97   8.00  20.00 1.00     1970     2055
##    y0[12]  12.46  12.00 3.68 2.97   7.00  19.00 1.00     1782     1982
##    y0[13]  20.61  20.00 4.97 4.45  13.00  29.00 1.00     1761     1797
##    y0[14]  19.10  19.00 4.70 4.45  12.00  27.00 1.00     2200     1998
##    y0[15]  14.14  14.00 3.88 4.45   8.00  21.00 1.00     1957     1752
##    y0[16]  12.80  13.00 3.64 2.97   7.00  19.00 1.00     1796     2006
##    y0[17]  10.89  11.00 3.49 2.97   5.00  17.00 1.00     1427     1653
##    y0[18]  12.52  12.00 3.73 2.97   7.00  19.00 1.00     1704     1960
##    y0[19]  17.50  17.00 4.43 4.45  11.00  25.00 1.00     2163     1923
##    y0[20]  10.57  10.00 3.25 2.97   5.00  16.00 1.00     1385     1821

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

ex13-1.stan

generalized linear mixed model, poisson regression

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
  real<lower=0> sr0;
  vector[N] r0;
}
model{
  r0~normal(0,sr0);
  for(i in 1:N)
    y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+r0[i]+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-1.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -30405.3 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99       10069.1     0.0268383       19129.8      0.2113     0.02113      145    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199       10309.9     0.0331786   8.55093e+13      0.3351      0.3351      342    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      206         10310    0.00319902   1.03472e+12       0.681       0.681      349    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__   10310.00
##    b0         2.84
##    b1         0.00
##    sr0        0.00
##    r0[1]      0.00
##    r0[2]      0.00
##    r0[3]      0.00
##    r0[4]      0.00
##    r0[5]      0.00
##    r0[6]      0.00
##    r0[7]      0.00
##    r0[8]      0.00
##    r0[9]      0.00
##    r0[10]     0.00
##    r0[11]     0.00
##    r0[12]     0.00
##    r0[13]     0.00
##    r0[14]     0.00
##    r0[15]     0.00
##    r0[16]     0.00
##    r0[17]     0.00
##    r0[18]     0.00
##    r0[19]     0.00
##    r0[20]     0.00
##    m0[1]      2.84
##    m0[2]      2.84
##    m0[3]      2.84
##    m0[4]      2.84
##    m0[5]      2.84
##    m0[6]      2.84
##    m0[7]      2.84
##    m0[8]      2.84
##    m0[9]      2.84
##    m0[10]     2.84
##    m0[11]     2.84
##    m0[12]     2.84
##    m0[13]     2.84
##    m0[14]     2.84
##    m0[15]     2.84
##    m0[16]     2.84
##    m0[17]     2.84
##    m0[18]     2.84
##    m0[19]     2.84
##    m0[20]     2.84
##    y0[1]     12.00
##    y0[2]     15.00
##    y0[3]     21.00
##    y0[4]     12.00
##    y0[5]     15.00
##    y0[6]     26.00
##    y0[7]     17.00
##    y0[8]     14.00
##    y0[9]      9.00
##    y0[10]    17.00
##    y0[11]    14.00
##    y0[12]    19.00
##    y0[13]    21.00
##    y0[14]    15.00
##    y0[15]    14.00
##    y0[16]    21.00
##    y0[17]    22.00
##    y0[18]    17.00
##    y0[19]    13.00
##    y0[20]     9.00
mcmc=goMCMC(mdl,data,wrm=500)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.8 seconds.
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 finished in 0.8 seconds.
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.9 seconds.
## Chain 4 finished in 0.8 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.8 seconds.
## Total execution time: 1.0 seconds.
seeMCMC(mcmc,'[m,r]')
##  variable    mean  median    sd   mad      q5      q95 rhat ess_bulk ess_tail
##    lp__   9984.76 9983.98 13.47 13.68 9963.32 10008.50 1.02       76       92
##    b0        2.23    2.23  0.03  0.03    2.19     2.28 1.00      702      722
##    b1        0.09    0.09  0.00  0.00    0.08     0.09 1.00      784      743
##    sr0       0.01    0.01  0.01  0.01    0.00     0.03 1.02       76       91
##    r0[1]     0.00    0.00  0.02  0.01   -0.03     0.03 1.01     4126      787
##    r0[2]     0.00    0.00  0.01  0.01   -0.03     0.02 1.01     2956      887
##    r0[3]     0.00    0.00  0.02  0.01   -0.03     0.03 1.01     2947      724
##    r0[4]     0.00    0.00  0.01  0.01   -0.03     0.02 1.00     3740      873
##    r0[5]     0.00    0.00  0.02  0.01   -0.02     0.03 1.01     2896      714
##    r0[6]     0.00    0.00  0.01  0.01   -0.02     0.02 1.01     3065      713
##    r0[7]     0.00    0.00  0.02  0.01   -0.03     0.03 1.01     2504      735
##    r0[8]     0.00    0.00  0.02  0.01   -0.02     0.02 1.00     3240      907
##    r0[9]     0.00    0.00  0.01  0.01   -0.02     0.02 1.00     2565      955
##    r0[10]    0.00    0.00  0.02  0.01   -0.02     0.02 1.00     2608      743
##    r0[11]    0.00    0.00  0.02  0.01   -0.02     0.03 1.01     3916      828
##    r0[12]    0.00    0.00  0.02  0.01   -0.02     0.02 1.01     3169      605
##    r0[13]    0.00    0.00  0.02  0.01   -0.03     0.03 1.01     3166      778
##    r0[14]    0.00    0.00  0.01  0.01   -0.02     0.02 1.01     3116      965
##    r0[15]    0.00    0.00  0.02  0.01   -0.02     0.02 1.00     2852      909
##    r0[16]    0.00    0.00  0.01  0.01   -0.02     0.02 1.00     3048      774
##    r0[17]    0.00    0.00  0.01  0.01   -0.02     0.02 1.00     3239     1041
##    r0[18]    0.00    0.00  0.02  0.01   -0.03     0.03 1.00     2791      694
##    r0[19]    0.00    0.00  0.01  0.01   -0.03     0.02 1.01     3670      862
##    r0[20]    0.00    0.00  0.01  0.01   -0.02     0.02 1.01     3273      959
##    m0[1]     2.24    2.24  0.03  0.03    2.19     2.29 1.00      842      852
##    m0[2]     2.25    2.25  0.03  0.03    2.20     2.30 1.00      844      998
##    m0[3]     2.45    2.45  0.02  0.02    2.41     2.48 1.00      986      943
##    m0[4]     2.30    2.30  0.03  0.03    2.25     2.35 1.01      881     1053
##    m0[5]     3.01    3.01  0.02  0.03    2.97     3.05 1.00     1774     1533
##    m0[6]     2.82    2.82  0.02  0.02    2.79     2.86 1.00     2348     1663
##    m0[7]     2.52    2.52  0.02  0.02    2.48     2.55 1.00     1198      968
##    m0[8]     2.94    2.94  0.02  0.02    2.90     2.98 1.00     2075     1350
##    m0[9]     3.05    3.05  0.03  0.03    3.00     3.09 1.00     1509     1359
##    m0[10]    2.76    2.76  0.02  0.02    2.73     2.79 1.00     2515      986
##    m0[11]    2.62    2.62  0.02  0.02    2.59     2.65 1.01     1927      901
##    m0[12]    2.52    2.52  0.02  0.02    2.48     2.55 1.00     1137     1189
##    m0[13]    3.02    3.02  0.03  0.02    2.97     3.06 1.00     1898     1388
##    m0[14]    2.95    2.95  0.02  0.02    2.91     2.99 1.00     1773     1269
##    m0[15]    2.66    2.66  0.02  0.02    2.63     2.69 1.00     1972     1231
##    m0[16]    2.56    2.56  0.02  0.02    2.52     2.59 1.00     1342     1128
##    m0[17]    2.39    2.38  0.03  0.02    2.35     2.43 1.00      922     1036
##    m0[18]    2.54    2.54  0.02  0.02    2.51     2.58 1.00     1111     1130
##    m0[19]    2.86    2.86  0.02  0.02    2.82     2.89 1.00     2567     1308
##    m0[20]    2.35    2.35  0.03  0.03    2.31     2.39 1.00      961      987
##    y0[1]     9.44    9.00  2.99  2.97    5.00    15.00 1.00     2104     2052
##    y0[2]     9.53    9.00  3.04  2.97    5.00    15.00 1.00     1947     2018
##    y0[3]    11.50   11.00  3.40  2.97    6.00    17.00 1.00     1989     1787
##    y0[4]     9.98   10.00  3.22  2.97    5.00    16.00 1.00     1898     2005
##    y0[5]    20.23   20.00  4.62  4.45   13.00    28.00 1.00     2159     1847
##    y0[6]    16.87   17.00  4.18  4.45   10.00    24.00 1.00     1887     1877
##    y0[7]    12.36   12.00  3.65  2.97    7.00    19.00 1.00     2082     2040
##    y0[8]    18.91   19.00  4.38  4.45   12.00    26.00 1.00     2031     2028
##    y0[9]    21.25   21.00  4.64  4.45   14.00    29.00 1.00     1919     1885
##    y0[10]   15.76   16.00  3.93  4.45   10.00    22.00 1.00     2012     2037
##    y0[11]   13.79   14.00  3.71  4.45    8.00    20.00 1.00     2045     2075
##    y0[12]   12.24   12.00  3.48  2.97    7.00    18.00 1.00     2077     2035
##    y0[13]   20.44   20.00  4.56  4.45   13.00    28.00 1.00     1987     1818
##    y0[14]   19.15   19.00  4.41  4.45   12.00    27.00 1.00     1997     2056
##    y0[15]   14.21   14.00  3.73  2.97    8.00    21.00 1.00     2104     1893
##    y0[16]   12.85   13.00  3.57  2.97    7.00    19.00 1.00     1981     1923
##    y0[17]   10.80   11.00  3.34  2.97    6.00    17.00 1.00     1989     1937
##    y0[18]   12.69   13.00  3.62  4.45    7.00    19.00 1.00     2127     1968
##    y0[19]   17.40   17.00  4.21  4.45   11.00    25.00 1.00     2019     1915
##    y0[20]   10.43   10.00  3.32  2.97    5.00    16.00 1.00     1962     1853

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')